function sac2()

sf = getmainselection;

if(sf > 0)
    pflag = getplotflag;
    QueMessage('SAC analysis, dual windows', 1); % clear the que
    for i = 1:length(sf)
        sac2_go(sf(i), pflag);
    end;
end;

function sac2_go(sf, pflag)

global DFILE ALLCH VOLTAGE CONTROL



isrc = 1;
t0 = 0;
tmax = 1000;
spike_thresh = 0;

X=[];
if(~isempty(CONTROL(sf).spike))
    X=getstructs(sf, 'spike.latency');
    X=X{:};
    ns = 0;
    for i = 1:length(X)
        ns = ns + length(X{i});
    end;
    fprintf(1, 'spike data for %s contains: %d trials with %d spikes\n', ...
        CONTROL(sf).filename, length(X), ns);
end;
if(isempty(X) | isempty(CONTROL(sf).spike)) % only read spikes if we really need to - otherwise the data is in the database already
    fprintf(1, 'Spike data not available: Reading %s\n', CONTROL(sf).filename);
    [DFILE, DPAR, err] = ana_setup(DFILE, sf);
    if(err ~= 0)
        return;
    end;

    if(~isempty(ALLCH))
        VOLTAGE = ALLCH{isrc};
    end;
    %[protocol, rate, records, pts, frec, lrec, time, TM, ZT, TL, VL] = analysis_setup2(DFILE, sf);
    [first_spike, first_isi, nr_spikes, splat]=find_spikes2(DFILE, VOLTAGE, t0, tmax, spike_thresh);
    X = {splat.latency};
    spike.latency = {splat.latency};
    spike.source = {splat.source};
    spike.fsl = first_spike;
    CONTROL(sf).spike=spike;
end;

% [DFILE, DPAR, err] = ana_setup(DFILE, sf);
% if(err ~= 0)
%     return;
% end;
% 
% if(~isempty(ALLCH))
%     VOLTAGE = ALLCH{isrc};
% end;
% %[protocol, rate, records, pts, frec, lrec, time, TM, ZT, TL, VL] = analysis_setup2(DFILE, sf);
% [first_spike, first_isi, nr_spikes, splat]=find_spikes2(DFILE, VOLTAGE, t0, tmax, spike_thresh);
% 
% 
% X = {splat.latency};
twin = 100;
binw = 0.025;
start1 = 100;
start2 = 550;
dur1 = 400;
dur2 = 400;

[y, yh1, hx1, mr1] = sac(X, twin, binw, start1, dur1);
[y, yh2, hx2, mr2] = sac(X, twin, binw, start2, dur2);
% perform gaussian fits on the histograms to identify the peaks and measure
% them. Using Molitor's routines (mrqfit).
% [...] = MRQFIT(F, P, X, Y, SIG, VP, LB, UB, IMAX, TOL)
% F = 'gaussian', P = [A0 A1 M1 S1 A2 M2 S2 ... AN MN SN]   
% SIG is sigma Y (def = 1); VP is vary array per parameter [0 or 1];
% LB, UB are upper and lower bounds. 
%
hx1c = hx1+0.5*(hx1(2)-hx1(1));
hx2c = hx2+0.5*(hx2(2)-hx2(1));
gpar1 = [1, 5, 0, 1]; % single gaussian centered on 0
scf = 1;
gulim1 = [0, 100, 10*scf, 20*scf];
gllim1 = [0, 0, 0.0, 0.01];
gvar1 =  [0 1 0 1];
nitermax = 100;
gpar2 = [gpar1 5 10*scf 1*scf];
gulim2 = [gulim1 100 50*scf 20*scf];
gllim2 = [gllim1 0 2*scf 0.1*scf];
gvar2 = [gvar1 1 1 1];

[fp11, chisq11, niter11, fitc11, err11, dep11] = mrqfit('gaussian', gpar1, hx1c*scf, yh1, [], gvar1, gllim1, gulim1, nitermax, []);
[fp21, chisq21, niter21, fitc21, err21, dep21] = mrqfit('gaussian', gpar2, hx1c*scf, yh1, [], gvar2, gllim2, gulim2, nitermax, []);

[fp12, chisq12, niter12, fitc12, err12, dep12] = mrqfit('gaussian', gpar1, hx2c*scf, yh2, [], gvar1, gllim1, gulim1, nitermax, []);
[fp22, chisq22, niter22, fitc22, err22, dep22] = mrqfit('gaussian', gpar2, hx2c*scf, yh2, [], gvar2, gllim2, gulim2, nitermax, []);


yg1 = gaussfunc(hx1c, fp11);
yg2 = gaussfunc(hx2c, fp12);
fwhmfac = 2*sqrt(2*log(2)); % note - log is ln.
fwhm1 = fwhmfac*fp11(4); % full width at half maximal height.
fwhm2 = fwhmfac*fp12(4);


hf = findobj('tag', 'SAC2');
if(isempty(hf))
    hf = figure;
    set(hf, 'Tag', 'SAC2');
end;
% generate the results structure: SAC
SAC.hx1 = hx1; % save the histograms.
SAC.hy1 = yh1;
SAC.hx2 = hx2;
SAC.hy2 = yh2;
SAC.mr1 = mr1;
SAC.mr2 = mr2;

% gaussian fit results:
SAC.Gfit1 = fp11; % parameters
SAC.niter1 = niter11; % iterations
SAC.chisq1 = chisq11; % fit error
SAC.err1 = err11; % parameter estimate error
SAC.dep1 = dep11; % dependency between parameters
SAC.fwhm1 = fwhm1;
SAC.NPH1 = max(yg1);
%
SAC.Gfit2 = fp12;
SAC.niter2 = niter12;
SAC.chisq2 = chisq12;
SAC.err2 = err12;
SAC.dep2 = dep12;
SAC.fwhm2 = fwhm2;
SAC.NPH2 = max(yg2);

CONTROL(sf).SAC = SAC; % save in the database.

figure(hf);
clf;


% text area
subplot('position', [0.1, 0.90, 0.8, 0.095]);
axis([0,1,0,1])
axis('off')
ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]',CONTROL(sf).filename, ...
    CONTROL(sf).recbeg, CONTROL(sf).recend, CONTROL(sf).protocol, date), 'Fontsize', 10);
set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
   text(0,0.6,sprintf('Solution:%-25s  ', CONTROL(sf).solution), 'FontSize', 8);
   text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
      CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);
   text(0,0.200,sprintf('Window1: %.1f-%.1f  mean rate: %.2f s/s  Window 2: %.1f-%.1f ms  mean rate: %.2f s/s', start1, start1+dur1, mr1, start2, start2+dur2, mr2), 'FontSize', 8);
   text(0, 0.000, sprintf('G1: A=%.2f (%.2f) S=%.3f (%.3f)  NPH = %.2f  FWHM1 = %.3f', ...
       fp11(2), err11(2), fp11(4), err11(4), SAC.NPH1, SAC.fwhm1), 'Fontsize', 8);
   text(0, -0.200, sprintf('G2: A=%.2f (%.2f) S=%.3f (%.3f) NPH2 = %.2f  FWHM2 = %.3f', ...
       fp12(2), err12(2), fp12(4), err12(4), SAC.NPH2, SAC.fwhm2), 'Fontsize', 8);
   

subplot('position', [0.1, 0.075, 0.8, 0.320]);
%bar(sqrt(hx1), yh1, 'histc');
semilogx(hx1c, yh1, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
hold on
semilogx(hx1c, yg1, 'r');
%semilogx(hx1c, fitc11, 'g');

u1 = get(gca, 'Ylim');
h1 = gca;
set(gca, 'Xlim', [0 100]);
xlabel('Delay (ms)');
ylabel('Normalized Peak Height');


subplot('position', [0.1, 0.480, 0.8, 0.320]);

%bar(sqrt(hx2), yh2, 'histc');
semilogx(hx2c, yh2, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
hold on
semilogx(hx2c, yg2, 'r');
%semilogx(hx2c, fitc12, 'g');

u2=get(gca, 'YLim');
h2 = gca;
set(gca, 'Xlim', [0 100]);
xlabel('Delay (ms)');
ylabel('Normalized Peak Height');
if(u1(2) > u2(2))
    u2(2) = u1(2);
end;
u2(1) = 0;
set([h1 h2], 'Ylim', u2);
y0 = u2(2);
x0 = 100;
axes(h1);
text(x0, y0, sprintf('Window 1 (%.1f - %.1f ms)', start1, start1+dur1), ...
    'horizontalalignment', 'right', 'verticalalignment', 'top', 'fontsize', 9);
box off
axes(h2);
text(x0, y0, sprintf('Window 2 (%.1f - %.1f ms)', start2, start2+dur2), ...
    'horizontalalignment', 'right', 'verticalalignment', 'top', 'fontsize', 9);
box off

if(pflag)
    orient landscape
print;
end;


function [y] = gaussfunc(x, fp)
%
% calculate a gaussian based on x and FP
%
y = fp(1) + (fp(2)/(fp(4)*sqrt(2*pi)))*exp(-((x-fp(3)).^2)/(2*fp(4)^2));

